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1.0  INTRODUCTION 


1.1  Background 

Th'S-'work  described  in  this  report  follows  that  performed  under  contract 
DACA39- 77-0-0004  and  presented  in  "Computer  Modelling  of  Jointed  Rock  Masses"  by 
Cundall  e£  al  (1978)*. 

In  the  previous  study  several  computer  programs  were  developed  as  exam¬ 
ples  of  the  "distinct  element  method",  which  is  an  explicit,  time-marching  pro¬ 
cedure  that  models  assemblies  of  discrete  blocks  or  particles  that  interact  mech¬ 
anically.  The  first  program  was  restricted  to  rigid  blocks,  and  was  a  translation 
into  Fortran  of  an  earlier  machine- language  code  that  used  interactive  graphics 
to  create  and  manipulate  blocks.  A  second  program  introduced  simple  deformability, 
whereby  each  block  has  three  degrees  of  freedom  to  deform  internally.  A  modified 
version  of  the  rigid  block  program  was  written  in  order  to  allow  blocks  to  split 
into  two,  using  a  simple  criterion  based  on  the  applied  loads  and  the  block  dim¬ 
ensions.  Finally,  an  experimental  program  was  written  in  which  blocks  could  be 
discretised  internally  into  finite-difference  triangles.  Such  blocks  are  termed 
"fully-deformable" . 

Apart  from  the  inconvenience  in  having  different  facilities  available 
in  different  programs,  the  computer  programs  were  written  with  no  particular 
emphasis  on  efficiency  or  flexible  data  structures,  since  the  intention  was  to 
demonstrate  and  test  some  new  formulations.  Furthermore,  it  has  become  apparent 
that  it  would  be  difficult  to  represent  certain  physical  behaviour,  such  as  fluid 
interaction,  edge-to-edge  contact  and  soft  comers,  without  major  overhaul  of 
the  program  logic  and  data  structures. 

These  considerations  prompted  the  program  development  described  in  this 

report . 

1.2  Scope  of  Present  Work 

A  completely  new  program,  UDEC+,  has  been  developed,  which  provides, 
in  one  package,  almost  all  of  the  capabilities  that  existed  separately  in  the 
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previous  programs.  The  principal  objective  was  to  write  a  program  that  would 
allow  rigid  blocks,  simply-deformable  blocks  and  fully-deformable  blocks  to  be 
mixed  together  in  one  numerical  simulation.  It  is  often  useful  to  mix  differ¬ 
ent  types  of  block  for  the  following  reasons:  near  a  free  surface  in  jointed 
rock,  the  movements  arise  predominantly  from  slip  and  opening  of  joints.  In 
these  regions  rigid  blocks  may  be  used,  for  maximum  efficiency  of  calculation. 

On  moving  away  from  a  free  surface,  towards  the  interior  of  a  rock  mass,  joint 
displacements  diminish  in  comparison  with  deformations  of  the  intact  rock,  and 
the  stress  distribution  is  determined  largely  by  the  elastic  properties  of  the 
rock.  Deformable  blocks  should  be  used  in  these  regions  to  represent  the  rock 
behaviour  correctly.  The  correct  boundary  conditions  are  especially  important 
for  dynamic  calculations  in  which  incident  waves  are  to  be  propagated  towards  a 
rock  structure,  and  reflected  waves  are  to  be  absorbed.  Again,  deformable  blocks 
are  needed  near  the  boundary  to  provide  the  correct  propagation  velocity  and  the 
correct  driving  impedence  for  absorbing  boundaries . 

Although  the  contractural  requirement  of  the  study  was  limited  to  the 
provision  of  a  program  incorporating  the  three  types  of  block  noted  above,  the 
opportunity  was  taken  to  re-examine  the  way  in  which  the  data  describing  the 
blocks  is  represented:  the  "data  structure".  The  data  structure  determines 
to  a  very  great  degree  how  easy  it  is  to  represent  diverse  physical  phenomena. 

Two  guiding  principles  influenced  the  choice  of  a  data  structure:  firstly  it 
was  assumed  that  the  best  data  structure  is  the  one  that  corresponds  in  a  topo¬ 
logical  sense  most  closely  to  the  physical  structure.  In  previous  programs, 
the  only  correspondence  was  between  stored  variables  and  physical  variables . 

The  new  program  actually  stores  the  variables  in  memory  in  a  way  that  has  iden¬ 
tical  topological  properties  to  the  physical  arrangement.  The  second  assumption 
that  was  kept  in  mind  during  the  design  of  program  UDEC  was  concerned  with  the 
way  in  which  computer  hardware  is  developing.  Typical  memory  sizes  have  in¬ 
creased  and  costs  have  decreased  by  orders  of  magnitude  over  the  last  decade; 
megabyte  memories  are  new  common  in  minicomputers,  and  are  just  becoming  avail¬ 
able  for  microcomputers .  Execution  speeds  have  also  increased,  but  at  a  slower 
rate.  Often,  when  writing  a  program,  it  is  possible  to  make  a  trade-off  between 
memory  use  and  execution  time:  for  example,  variables  may  be  saved  to  avoid  re¬ 
calculation  within  a  loop.  The  philosophy  of  reducing  execution  time  at  the 
expense  of  greater  memory  requirements  has  been  adopted  throughout. 
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The  numerical  formulations  of  the  various  block  types  are  almost 
identical  to  those  documented  in  the  previous  report  by  Cundall  et  al  (loc  cit) . 
The  formulations  are  not  repeated  in  this  report,  and  the  reader  is  encouraged 
to  read  the  previous  report  in  conjunction  with  this  one.  Any  differences, 
such  as  rounded  comers,  are  described  in  Chapter  2,  which  also  documents  other 
changes  and  additions  arising  from  the  new  data  structure.  The  data  structure 
itself  is  documented  in  detail  in  Chapter  3. 

Program  UDEC  allows  interaction  between  three  block  types:  rigid, 
simply-deformable  and  fully-deformable .  Of  the  six  possible  interactions 
(e.g.  rigid/rigid,  rigid/simply-deformable  etc.)  five  require  (due  to  their  basic 
formulation)  that  contacts  shall  consist  of  finite  stiffnesses.  Only  the  inter¬ 
action  between  fully-deformable  blocks  has  the  possibility  of  rigid  contact 
between  opposing  grid-points.  Rezoning  is  necessary  in  this  case  to  keep  con¬ 
tacting  grid-points  opposite  one  another.  It  was  decided  to  omit  re zoning  from 
program  UDEC,  and  impose  a  restriction  that  all  contacts  be  of  finite  stiffness. 

In  practice,  an  arbitrarily-large  stiffness  may  be  used,  but  at  the  expense  of 
a  very  small  time  step.  It  was  felt  that  the  complication  of  rezoning  was  in¬ 
appropriate  in  a  program  that  was  primarily  intended  to  model  rock  blocks  with 
finite-stiffness  joints,  and  that  its  inherent  inefficiency  could  not  be  justified 
by  the  very  small  number  of  cases  for  which  it  might  be  used. 


2.0  ADDITIONS  AND  CHANGES  TO  PREVIOUS  WORK 


There  are  a  number  of  problems  associated  with  existing  distinct  ele¬ 
ment  programs  for  modelling  assemblies  of  angular  blocks.  The  more  important 
of  these  are  as  follows: 

a.  The  system  of  "boxes"  used  for  coarse  classification  of 
blocks  is  quite  efficient  when  the  blocks  are  of  a  similar 
size.  For  blocks  that  are  much  larger  them  the  box 
dimension,  much  time  is  taken  searching  all  the  boxes 
overlapping  the  block  and  depositing  and  moving  entries 
among  these  boxes .  For  blocks  that  cure  small  compared 

to  the  box  size,  many  potential  contacts  must  be  tested 
within  each  box,  although  only  a  few  will  be  accepted. 

b.  It  is  necessary  to  perform  updates  (searches  for  new 
contacts)  globally  in  all  programs  except  RBMC,  owing  to 
the  difficulty  of  guaranteeing  that  a  given  corner  will 
always  locate  a  nearby  edge. 

c.  There  is  a  limitation  on  the  coordinate  system  that  can 
be  chosen,  so  that  the  dimensions  of  a  given  block  system 
must  be  scaled  before  running  a  problem.  This  restriction 
is  made  because  the  magnitudes  of  block  coordinates  are 
used  to  trigger  re-boxing  of  pointers  within  the  box  system. 

d.  Blocks  cannot  be  deleted,  except  in  the  original  machine- 
language  version  of  the  rigid-block  program. 

e.  "Hang-ups"  occur  when  two  blocks  overlap  by  an  arbitrarily- 
small  amount,  because  block  comers  are  assumed  to  be  sharp 
and  infinitely  strong. 

f.  Contact  is  always  between  an  edge  and  a  corner;  edge-to-edge 
contact  is  only  approximately  represented  by  two  edge-to-corner 
contacts . 
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g.  When  the  edge  of  one  block  slides  past  the  edge  of  another, 
as  shown  below,  there  is  a  discontinuity  in  contact  force 
because  one  edge-to-corner  contact  must  break  before  another 
comer- to-edge  contact  can  form.  Stored  energy  is  lost 
because  a  load-carrying  contact  is  suddenly  deleted. 


FIXED 


=i 

CORNER -T0- 
EDGE  CONTACT 


EDGE- TO-  I 
CORNER  CONTACT 


All  the  problems  noted  above  are  eliminated  in  the  new  formulation 
of  the  distinct  element  method  described  in  this  report.  The  main  changes  are 
to  the  data  structures,  and  in  the  fact  that  rounded  comers  are  used.  Several 
specific  innovations  are  described  below. 

2.1  Detection  of  Contacts 

Chapter  3  describes  the  new  data  structues  in  detail,  and  the  way 
in  which  their  topological  properties  correspond  closely  with  those  of  the 
physical  system  of  blocks  that  they  represent.  It  is  this  correspondence  that 
allows  the  old  "box"  scheme  to  be  dispensed  with:  the  connectivity  of  the  physical 
system  is  built  into  the  data  structure,  which  means  that  potential  contacts  with 
a  given  block  may  be  detected  by  a  local  examination  of  the  linked-list  network 
surrounding  the  block.  However,  for  the  scheme  to  work  well,  there  must  be  a 
well-developed  connectivity;  the  main  application  will  be  to  systems  of  blocks. 
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each  of  which  is  very  near  several  other  blocks.  This  is  likely  to  be  true 
for  program  UDEC,  which  has  the  main  applications  in  modelling  jointed  rock 
masses,  where  blocks  are  all  touching  initially. 


2.2  Update  Triggering 

An  "update"  is  defined  here  as  a  scan  of  some  region  to  determine  if 
new  contacts  should  be  made  or  to  erase  unwanted  contacts.  As  illustrated  be¬ 
low,  a  new  contact  can  only  arise  physically  within  a  closed  region  between 
blocks,  called  a  "domain". 


For  this  reason,  an  update  is  limited  to  an  individual  domain,  and  is 
triggered  by  significant  relative  motion  occurring  within  the  domain.  A  scheme 
has  been  adopted  whereby  a  fictitious  displacement  is  accumulated  for  each  domain. 
This  displacement  is  related  to  the  relative  motion  that  has  taken  place  in  the 
domain  since  the  previous  update,  and  is  used  to  trigger  the  next  update  when 
the  displacement  exceeds  a  certain  tolerance.  At  each  time-step,  the  greatest 
relative  x-  and  y-velocities  are  recorded  between  any  two  corners  within  the  do¬ 
main  (including  corners  that  are  part  of  a  contact) .  The  fictitious  displace¬ 
ment  is  then  accumulated  as  follows: 


where  um  and  um  are  the  greater  relative  velocities,  u^  is  the 
x  y 

fictitious  displacement  and  At  is  the  time-step.  An  update 
is  done  for  the  domain  when 

uf  >  35. 

/2 

where  T  is  the  tolerance  on  making  contacts  ("contacts"  are  established 
m 

before  the  two  blocks  actually  touch:  T  is  the  distance  between 

m 

potentially  contacting  points) .  The  criterion  for  triggering 
updates  given  above  ensures  that  contacts  are  always  detected 
before  physical  contact  is  made. 

2.3  Rounded  Corners 

In  program  UDEC,  blocks  behave  mechanically  as  if  each  corner  consists 
of  an  arc  of  a  circle;  the  arc  is  tangent  to  the  two  adjoining  edges.  The 
circle  is  defined  by  specifying  the  distance  from  the  block  comer  to  the  inter¬ 
section  points  of  the  circular  arc  with  the  adjoining  edges.  This  procedure 
seems  physically  more  reasonable  than  a  specification  of  a  constant  radius  for 
all  corners,  since  sharp  comers  would  be  considerably  truncated  if  the  same 
radius  was  used  for  all  comers . 


CONSTANT  RADIUS 


CONSTANT  DISTANCE  TO  CORNER 


The  distance,  D,  may  be  specified  by  the  user,  and  applies  to  all 
blocks.  It  is  stored  as  Fortran  variable  OTOL.  Radius  and  circle  centre 
are  calculated  for  each  corner  as  follows: 


R 


X  . 

CORNER  1 


ZL  -  {(x£  -  x.)  <*J  -  X.)}*2 

„R  i,  R  ,  ,  R 
Z  =  <(x,  -  x. )  (x.  -  x. )  > 
l  i  i  i  i  f 

unit  vectors: 

U^  =  (xj  -  xi)/ZL 

R  ,  R  i  j„R 
u  =  (x±  -  X  )/Z 

tangent  unit  vectors: 


by  vector  addition 


DuL  +  Rtf  +  Rtf  =  Duf 
i  i  i  i 


therefore 


R 


or 


R 


R 


(1) 


P(u2  ~  u2} 
uj  +  uj  (2) 


Either  form  (1)  or  form  (2)  may  be  used,  depending  on  the  magnitude 
Li  R  Li  R 

of  u^  +  u2  compared  to  u^  +  u^.  The  expression  with  the  largest  denominator 
is  used  in  UDEC. 

The  circle  centre  is  found  from: 

C  L  L 

Xi  =  Xi  +  DUi  +  Rti 

In  the  present  form  of  UDEC,  corners  are  assumed  to  be  rounded  only 
for  the  purpose  of  calculating  contact  mechanisms;  block  mass  and  moment  of 
inertia  are  calculated  on  the  assumption  that  comers  are  angular.  This  assump¬ 
tion  is  not  essential,  and  the  program  can  be  changed  quite  easily. 

2.4  Edge-to-edge  contact 

Chapter  3  contains  a  description  of  the  data  structures,  and  in  partic¬ 
ular  the  linked- list  representation  of  a  domain.  A  domain  is  a  closed  area 
bounded  by  blocks  and  contacts:  the  associated  linked- list  is  a  circular  chain 
of  contacts  and  comers  encountered  during  an  anti-clockwise  scan  of  the  domain 
boundary . 


Only  two  types  of  contact  are  needed  by  the  data  structure  for  represent¬ 
ing  a  system  of  blocks:  comer- to-corner  contacts  and  edge-to-comer  contacts. 


These  will  be  termed  "numerical  contacts".  Physically,  however,  edge-to-edge 
contact  is  important  because  it  corresponds  to  the  case  of  a  rock  joint  closed 
along  its  whole  length.  A  physical  edge-to-edge  contact  corresponds  to  a  domain 
with  exactly  two  numerical  contacts  in  its  linked-list.  When  an  edge-to-edge 
contact  is  recognised  in  this  way  during  a  domain  scan,  it  is  treated  differently 
from  other  contacts  in  respect  of  its  physical  behaviour.  For  example ,  it  is 
more  appropriate  to  express  the  constitutive  model  for  such  a  contact  in  terms 
of  stresses  rather  them  forces.  Stresses  cem  be  evaluated  by  the  program  since 
the  length  of  the  physical  contact  is  known  from  the  distance  between  the  two 
numerical  contacts. 

Memy  types  of  constitutive  model  for  edge-to-edge  contact  may  be  con¬ 
templated.  The  program  provides  the  displacements  at  either  end  of  the  joint, 
and  the  model  must  furnish  the  average  normal  and  shear  stress,  and  the  line  of 
action  of  the  resultant  forces.  The  simplest  consitutive  model  is  as  follows: 


_ Li  B 

v  y  Un 
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Aa  =  f  (a  ,  Au  ,  Au  ) 
n  n  n  n  s 


where:  Au  ■  —  (AuA  +  AuB) 
n  2  n  n 


At  =  f  (a  ,  t  ,  Au  ,  Au  ) 
s  n  ns 


Au  =*  (AuA  +  AuB) 
s  2  s  s 


a  =  a  +  Aa 
n  n  n 


A  A 

Au  ,  Au  =  normal  and  shear  displacement 
n  increments  at  end  A  of  joint; 


t  =  t  +  At 


B  B 

Au  ,  Au  =  same,  at  end  B  of  joint 
n  s 


a  =  normal  stress 
n 


t  =  shear  stress 


prime  (  )  denotes  new  value. 


f  ,  f  =  non  linear  functions 
n  s 


F  =»  p.a 
n  n 


p  =  joint  length 


F  =  p.a 
s  s 


F  ,F  =  overall  normal  and  shear  forces 
n  s 


FA  *  r  F 
n  n 


FA  =  r  F 
s  s 


FB  =  (l-r)F  FB  =  (l-r)F 
n  ns  s 


n 


A  .  B 
u  +  u 
n  n 


In  the  procedure  presented  above,  the  line  of  action  of  the  resultant 
forces  is  determined  by  making  the  ratio  of  end  forces  equal  to  the  ratio  of  end 
displacements : 
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Any  other  scheme  may  be  substituted:  for  example,  the  angle  between 
the  two  block  faces  may  be  used  to  determine  the  line  of  action,  or  the  centroid 
of  the  overlap  area  may  be  used.  In  both  cases,  the  program  already  provides 
the  necessary  geometrical  data. 

2.5  Contact  re-classification 

One  advantage  of  rounded  comers  is  that  there  are  no  abrupt  transitions 
as  blocks  in  contact  slide  past  one  smother.  In  the  sequence  shown  below,  the 
same  contact  is  re-classified  successively  as  corner/edge,  comer/ comer,  edge/ 
corner . 


In  previous  block  programs,  the  same  sequence  would  involve  the  initial 
contact  breaking,  and  a  new  one  being  created.  The  data  structures  of  UDEC  allow 
the  same  contact  to  be  retained  through  the  whole  sequence.  The  only  change  is 
that  a  code  number  denoting  the  classification  of  the  contact  is  updated  as  the 
blocks  move. 

2.6  Calculation  sequence 

In  all  explicit,  time-marching  schemes  the  main  calculation  cycle  con¬ 
sists  of  applying  the  law  of  motion  to  all  mass-points  followed  by  the  calculation 
of  force  increments  from  displacement  increments  for  all  spring-like  elements 
(contacts,  continuum  zones,  fluid  cells,  etc.).  Program  UDEC  follows  this  gen¬ 
eral  scheme  although  it  is  complicated  by  the  need  to  make  and  break  contacts, 
and  the  fact  that  many  different  types  of  element  are  allowed  to  interact. 
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For  all  blocks : 


.  accelerate  centroids  from  force-sums. 

.  calculate  strain-rates  for  simply-deformable  blocks  from  applied 
stresses . 

.  accelerate  grid-point  masses  from  internal  and  boundary  forces 
for  fully-deformable  blocks. 

.  update  comer  velocities  and  displacements . 

.  apply  new  relative  velocities  to  surrounding  contacts. 

.  reset  force-sums  to  zero. 

For  all  domains; 

.  accumulate  fictitious  domain  displacements  and  update  domain  if 
necessary  (i.e.  make  and  break  contacts) . 

.  compute  incremental  pore  pressures  from  velocities  around  domain; 
apply  resulting  forces  to  blocks. 

For  all  contacts: 

.  update  contact  forces  from  known  relative  contact  velocities  using 
constitutive  model. 

.  accumulate  centroid  force-sums,  grid-point  force-sums  and  applied 
stresses  for  simply-deformable  blocks  from  contact  forces. 

.  allow  fluid  flow  between  the  2  domains  on  either  side  of  contact; 
update  pore  pressures. 

For  all  zones; 

.  compute  strain  rates;  hence  new  stresses;  hence  grid-point  forces. 

2.7  Creation  of  blocks  and  joints  by  splitting 

The  main  application  of  program  UDEC  will  be  to  model  in-si tu  rock  masses. 

It  is  convenient  to  generate  the  numerical  system  of  rock  blocks  by  specifying 
the  joints  (discontinuities)  rather  than  individual  blocks.  Not  only  is  the 
input  data  less  voluminous,  but  the  joint  properties  can  be  specified  independently 
of  the  rock  properties.  Consequently  the  principal  way  of  creating  a  system  of 
rock  blocks  in  UDEC  is  by  splitting  existing  blocks,  and  deleting  blocks  that  are 
not  needed.  Almost  any  geometrical  arrangement  may  be  created  in  this  way. 

The  subroutine  that  splits  blocks  may  also  be  called  dynamically,  while  UDEC  is 
running,  so  that  blocks  may  fracture  as  a  function  of  the  load  applied  to  them. 


2.8  Pore-pressure  calculations 

The  framework  for  representing  pore-pressure  generation  and  dissipation 
exists  in  UDEC  but  no  tests  of  the  capability  have  been  made  at  the  time  of 
writing.  The  data  structure  provides  a  network  of  inter-connected  domains  and 
contacts  that  can  be  regarded  as  reservoirs  and  pipes  respectively.  For  each 
domain,  increments  in  fluid  pressure  can  be  calculated  from  the  known  incremental 
displacements  of  the  block  comers  that  constitute  the  boundary  to  the  domain . 

The  block  forces  arising  from  the  fluid  pressures  can  be  calculated,  since  the 
coordinates  of  each  point  around  a  domain  are  known.  During  a  contact  scan, 
the  pressures  on  either  side  of  each  contact  can  be  relaxed,  by  assuming  that 
the  constriction  corresponding  to  the  contact  has  a  certain  hydraulic  conduc¬ 
tivity.  In  this  way  pore-pressures  may  be  continuously  generated  and  dissipated 
as  the  blocks  move. 

2 .9  Fully-deformable  blocks 

As  discussed  in  the  Introduction,  no  re-zoning  is  performed  when  fully- 
deformable  blocks  slide  over  one  another.  The  use  of  finite  contact  stiffnesses 
between  blocks  enables  the  calculations  for  adjacent  blocks  to  be  decoupled; 
they  communicate  through  a  common  boundary  force.  The  mass-points  (also  called 
grid-points)  on  the  boundary  of  a  fully-deformable  block  are  accelerated  in  the 
normal  way  during  each  time  step  in  proportion  to  the  sum  of  the  forces  arising 
both  from  within  the  block  and  from  any  contacts  that  exist  at  the  points.  Rel¬ 
ative  contact  velocities  are  updated  from  the  velocities  of  the  block  boundaries 
on  either  side  of  each  contact. 

The  general  scheme  outlined  above  is  complicated  by  the  fact  that  grid- 
points  from  opposing  blocks  need  not  coincide;  edge-to-coraer  contact  is  allowed. 

In  order  to  calculate  forces  and  displacements  at  boundary  locations  other  than 
grid-points,  it  is  assumed  that  each  boundary  segment  between  grid-points  acts  as 
a  rigid  bar  with  prescribed  end  velocities:  the  velocities  are  assumed  to  vary 
linearly  along  the  bar.  This  assumption  is  consistent  with  the  assumption  that 
the  triangular  zones  are  of  the  constant-strain  type .  Any  forces  acting  on  the 
tar  are  distributed  to  the  two  ends,  while  maintaining  moment  and  force  equilibrium. 


14 


grid-point 


.B 

u. 

i 


applied 

force 


uA,  uB  =  end  velocities 

i  l 


f 

I 

FULLY -DEFORMABLE  VELOCITY  ASSUMPTION  FORCE  ASSUMPTION 

’  BLOCK 

A  further  complication  is  caused  by  the  rounded  corners .  In  the  case 
of  rigid  and  simply-deformable  blocks,  forces  and  velocities  are  transmitted 
exactly  at  the  point  of  contact.  For  fully-deformable  blocks,  the  calculation 
to  do  this  becomes  lengthy,  and  in  some  cases  ambiguous  when  several  zones  meet 
at  the  same  grid-point.  The  present  version  of  UDEC  contains  an  approximate, 
but  rapid  calculation  when  transmitting  forces  and  velocities  at  comers  that 
form  part  of  a  contact:  the  contact  force  is  assumed  to  act  at  the  corner  of 
the  block,  and  not  at  the  true  contact  point  found  by  assuming  rounded  corners, 
i  Similarly  the  velocity  transmitted  to  the  contact  is  taken  as  the  corner  velocity 

|  even  though  contact  is  located  somewhat  inside  the  corner-point.  The  rounded 

comers  still  operate  as  intended,  since  the  contact  normals  and  the  relative 
block  positions  for  contact  are  correct,  but  the  slight  error  in  location  of 
contact  forces  and  velocities  introduces  a  moment  and  rotation  error,  respectively. 

J  It  will  be  necessary  to  gain  experience  using  the  program  to  find  out  whether  the 

errors  are  significant  in  practice:  if  they  are  too  large,  the  program  cam  be 
modified,  but  at  the  expense  of  efficiency. 


3 .0  DATA  STRUCTURES 


3.1  Main  Structure 

There  are  many  subsidiary  pointers  and  lists  that  help  in  retrieving 
data  quickly,  but  the  main  data  structure  is  described  as  follows.  Each  block 
has  a  circular  linked- list  that  corresponds  to  its  boundary.  Each  corner  on 
the  boundary  is  represented  by  am  array  of  words  that  stores  the  coordinates  of 
the  corner,  the  velocity  and  several  other  items.  The  corner  arrays  are  linked 
together  consecutively  in  the  clockwise  direction.  When  two  blocks  come  into 
contact  their  corner  lists  are  broken  at  the  points  of  contact  and  an  array  of 
words  (corresponding  to  a  contact)  is  inserted  into  the  break,  such  that  the  con¬ 
tact  array  is  common  to  the  comer  lists  of  both  blocks.  The  process  of  making 
a  contact  is  illustrated  in  Pig.  1.  It  should  be  noted  that  in  all  diagrams, 
computer  variables  and  symbols,  the  following  letters  are  used  in  referring  to 
the  various  entities: 

C  contact  D  domain 

P  comer  Z  zone 

B  block  G  grid-point 

The  small  boxes  labelled  "comer  array"  and  "contact  array"  are  placed 
on  the  diagram  near  the  physical  locations  that  they  represent,  but  in  the  program 
UEEC  they  exist  as  groups  of  contiguous  memory  addresses  contained  in  the  single 
Fortran  array  IA (  ) .  The  actual  locations  of  the  groups  in  IA  (  )  are  unimportant 
and  arbitrary,  since  the  data  in  the  groups  are  retrieved  by  following  the  appro¬ 
priate  pointers. 

The  example  shown  in  Fig.  1  illustrates  that  there  are  two  paths  leading 
away  from  a  single  contact.  Depending  on  the  scanning  strategy  used,  it  is 
possible  to  follow  either  block  boundaries  or  the  spaces  (domains)  between  blocks. 
Fig.  2  illustrates  the  two  possibilities  for  the  same  data  structure. 

Fig.  2(a)  shows  an  anti-clockwise  scan  that  traces  the  boundary  of  an 
inter-block  domain;  Fig.  2(b)  shows  a  clockwise  scan  that  traces  a  block  boundary. 
The  first  type  of  scam  is  useful  when  calculating  pore  pressures  from  changes  in 
pore  volumes,  while  the  second  type  of  scam  is  used  when  up-dating  the  velocities 
of  comers  and  contacts  around  a  block  from  the  motion  of  the  block  itself. 
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FIG.  I  :  CHANGES  IN  LINKED  LISTS  THAT  OCCUR  WHEN  TWO 
SLOCKS  COME  INTO  CONTACT 
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a)  „  /  SCAN  AROUNO  THE  DOMAIN  BETWEEN  BLOCKS 


b)  jS  EXAMPLE  OF  SCAN  AROUNO  A  BLOCK 

FIG.  2  '  TWO  TYPES  OF  SCAN  THAT  ARE  POSSIBLE  WITH  THE 
SAME  OATA  STRUCTURE 
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The  data  structure  outlined  above  is  better  them  those  used  in  previous 
distinct  element  programs  for  the  following  reasons: 

a)  each  block  has  direct  access  to  all  its  contacts,  and  thereby  to  all 
neighbouring  blocks;  the  creation  and  deletion  of  contacts  is  made 
easy; 

b)  closed  domains  between  blocks  are  defined  without  additional  computing 
overhead;  pore-pressure  increments  may  be  calculated  from  known 
velocities  around  the  domain  boundaries ; 

c)  each  contact  has  direct  access  to  the  two  domains  on  either  side  of 
it;  hence  the  calculation  of  fluid  flow  between  domains  is  made  easy. 

d)  edge-to-edge  contact  may  be  detected  simply  by  noting  when  a  domain 
volume  collapses  to  zero. 

e)  The  search  for  new  contacts  can  be  triggered  locally,  by  monitoring 
relative  displacements  within  each  domain;  the  part  of  a  block  that 
intrudes  into  a  domain  can  only  touch  another  block  in  the  domain. 

The  observations  made  above  support  the  contention  made  in  the  Intro¬ 
duction  that  the  best  data  representation  is  that  which  corresponds  most  closely 
to  the  physical  structure:  any  changes  and  interactions  that  occur  physically 
can  be  treated  in  exactly  the  same  way  in  computer  memory. 

However,  the  data  structure  described  above  is  not  suitable  for  repre¬ 
senting  assemblies  of  particles  that  are  widely-separated,  or  particles  that  move 
with  high  velocities.  The  linked-list  scheme  assumes  that  the  connectivity  of 
the  system  changes  gradually,  and  that  the  domains  between  particles  are  clearly 
defined  by  a  closed  loop  of  contacts.  The  scheme  can  handle  isolated  cases 
where  individual  particles  lose  contact  with  their  neighbours,  by  using  "virtual 
contacts",  which  are  fictitious  links  between  particles.  But  the  computer  program 
is  not  designed  for  multiple  virtual  contacts;  in  this  case  it  becomes  inefficient. 

3.2  Data  arrays 

The  main  data  structure  described  above  involves  two  types  of  data  array: 
corner  and  contact  arrays.  Program  UDEC  uses  a  number  of  other  arrays  (or  groups 
of  contiguous  memory) .  The  complete  set  is  as  follows: 
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1 .  Block 

2 .  Domain 

3 .  Comer 

4 .  Contact 

5 .  SDEF  extension 

6 .  Zone 

7.  Grid-point 

8.  Domain  extension 

Appendix  III  lists  the  complete  contents  of  these  data  groups.  Groups 
3  and  4  have  already  been  described,*  Group  1  stores  data  for  each  block,  such  as 
centroid,  velocity,  material  type,  and  constitutive  type  number.  Group  2  stores 
the  pore-pressure  for  a  domain  and  an  accumulated  displacement  that  is  used  to 
trigger  a  contact  update  for  the  domain.  Group  5  contains  the  data  necessary 
for  simply-deformable  blocks,  such  as  stresses  and  s train- rates .  This  group  is 

linked  to  the  corresponding  block  data  (Group  1)  by  an  extension  pointer  that  is 
zero  in  the  case  of  a  rigid  block .  Groups  6  and  7  store  geometrical  and  connec¬ 
tion  data  for  fully-deformable  blocks ,  where  the  blocks  sure  discretised  into 
triangular  zones  and  grid-points.  Again,  the  linkage  from  the  corresponding  Group 
1  is  via  the  extension  pointer.  Group  8  stores  data  when  edge-to-edge  contacts 
are  detected.  Each  Group  8  data  block  is  linked  to  the  corresponding  domain  data 
(Group  2)  via  an  extension  pointer. 

3.3  Support  Structures 

A  number  of  other  lists  and  pointers  are  used,  so  that  data  can  be 
accessed  rapidly  and  conveniently,  particularly  when  cycling  through  the  main  cal¬ 
culation  loop. 

All  blocks,  domains,  contacts,  zones,  and  grid-points  are  linked  together 
into  five  separate  lists.  Pointers  to  the  starting  address  of  each  list  are  pro¬ 
vided;  Fig.  3  shows  this  arrangement  schematically. 

In  order  to  scan  through  block  comers  rapidly,  and  to  access,  con¬ 
veniently,  the  cornels  to  either  side  of  a  given  comer,  a  "reverse  list”  is 
provided  for  each  block  that  links  the  comers  of  a  block  together  in  the  anti¬ 
clockwise  direction.  It  will  be  recalled  that  the  normal  clockwise  list  of 
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FIG.  3  :  LINKED  LISTS  FOR  MAIN  DATA  ARRAYS 
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corners  could  be  broken  at  any  number  of  places  to  include  a  contact  array.  The 
reverse  list  has  no  such  interruptions:  it  simply  links  comers.  Fig.  4  illus¬ 
trates  the  reverse  list. 

Each  contact  joins  two  blocks  and  two  domains;  pointers  sure  provided 
in  the  contact  array  to  these  four  addresses.  Two  further  words  in  the  contact 
array  point  to  the  two  comer  lists  circulating  around  the  two  blocks.  It  is 
important  to  establish  an  ordering  convention  for  contact  pointers,  as  several 
exits  are  possible  after  entering  a  contact  during  a  scan.  Fig.  5  gives  this 
convention.  For  example,  when  entering  from  Block  1  (with  block  number  equal 
to  the  contents  of  the  word  with  offset  KB1) ,  the  domain  to  the  left  will  be  found 
from  the  word  with  offset  KDl. 

A  further  convention  is  observed  in  the  ordering  of  contacts  and  comer 
data  in  the  anti-clockwise  list  running  around  a  block  boundary.  All  contacts 
appear  in  the  list  in  the  order  in  which  they  occur  physically.  Since  comers 
are  rounded,  the  actual  location  of  a  comer  is  ambiguous,  particularly  if  there 
are  several  contacts  on  one  corner.  The  convention  is  adopted  whereby  the  comer 
data  array  is  always  placed  after  all  the  contact  data  arrays  for  that  comer. 

Fig.  7  illustrates  this  convention.  The  convention  is  convenient  because  a  given 
contact  can  always  locate  its  associated  comer (s)  by  following  the  linked  list 
until  a  comer  is  found. 

Several  data  arrays  contain  additional  pointers  not  noted  above:  a  com¬ 
plete  list  is  given  in  Appendix  III  for  each  data  type.  As  an  example,  each  group 
of  corner  data  for  a  fully-deformable  block  includes  a  pointer  to  the  data  group 
of  the  corresponding  grid-point.  Conversely,  the  grid-point  data  incorporates  a 
reverse  link  to  the  comer  data. 

3.4  Memory  Management 

All  requests  for  memory  are  handled  by  the  subroutine  FIND.  FIND  is 
given  the  number  of  words  required,  and  responds  by  returning  the  address  of  a 
newly  allocated  data  group,  if  memory  is  available.  Fresh  memory  is  used  to 
provide  space  for  a  new  array  unless  a  group  with  the  correct  number  of  words  has 
previously  been  returned.  This  reclamation  of  memory  is  made  possible  by  main¬ 
taining  a  linked  list  of  redundant  memory  groups.  The  list  is  constructed  by 
subroutine  LOSE  and  is  pointed  to  by  variable  JUNK,  and  its  structure  is  illus¬ 
trated  in  Fig.  6. 
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No  attempt  is  made  to  subdivide  memory  groups  when  smaller  numbers  of 
words  are  required,  but  this  technique  could  be  used  if  desired.  The  implemen¬ 
tation  of  a  suitable  scheme  should  not  be  too  difficult,  since  all  memory  manage¬ 
ment  is  handled  through  only  two  routines. 

3.5  Access  to  data  in  the  main  array 

Integers  and  real  numbers  are  both  stored  in  the  main  array,  IA(  ) . 

In  order  to  save  space  on  small  computers,  the  program  CJDEC  has  been  written  in 
such  a  way  that  the  word  length  for  integers  can  differ  from  the  word  length  for 
real  numbers.  The  following  schemes  are  possible,  for  example: 


Integers 

Reals 

16  bits 

32  bits 

32  bits 

32  bits 

16  bits 

64  bits 

32  bits 

64  bits 

64  bits 

64  bits 

This  flexibility  has  been  achieved  by  adopting  a  particular  convention 
when  accessing  data  from  the  main  array.  Integers  are  accessed  directly,  as 
follows: 

N  =  IA(IADDR) 

Real  numbers  are  always  accessed  via  a  subroutine  call: 

CALL  SUB (IA (IADDR) ) 

SUBROUTINE  SUB  (A) 

COMMON/  /  KOFF 

DIMENSION  A (1) 

R  -  A (KOFF) 

Here  R  is  the  required  real  .number  and  KOFF  is  the  offset  between  the 
calling  address  IADDR  and  the  location  of  the  number  R.  IADDR  is  typically  the 

pointer  to  the  start  of  a  particular  data  array,  and  KOFF  is  the  offset  corres¬ 
ponding  to  the  particular  variable  being  accessed. 
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As  set-up  at  present,  the  program  takes  the  space  of  two  Integers  to 
store  one  real  number.  This  can  be  changed  by  giving  different  values  to  the 
data  offsets,  which  are  variables  starting  with  the  letter  K.  Appendix  IV 
gives  further  details . 
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4.0  USE  OF  PROGRAM 


The  operation  of  program  UDEC  is  very  similar  to  that  of  the  previous 
programs  RBM  and  SEEM,  but  the  user  has,  in  general,  more  flexibility  in  the 
sequence  of  operations  that  may  be  requested.  Appendix  I  lists  the  input  com¬ 
mands  and  provides  a  description  of  each.  Appendix  II  contains  the  input  and 
output  from  some  complete  runs.  The  runs  illustrate  the  use  of  all  available 
commands ,  as  well  as  providing  sample  problems,  which  can  be  used  to  check 
versions  of  UDEC  set  up  on  other  computers. 

4.1  Creation  of  Blocks 

The  creation  of  the  assembly  of  blocks  differs  greatly  from  that  of 
previous  programs.  Instead  of  specifying  each  block  individually,  a  single, 
large  block  is  created  and  then  subdivided  repeatedly  into  many  small  blocks. 

Blocks  may  be  deleted  at  any  stage  in  order  to  generate  assemblies  of  complex 
geometry.  Blocks  may  be  split  and  deleted  even  after  cycling  has  started. 

4.2  Block  Types 

Blocks  created  initially  are,  by  default,  rigid.  Individual  blocks 
or  groups  of  blocks  subsequently  may  be  changed,  before  running,  to  simply- 
deformable  or  fully-deformable  blocks.  The  program  also  allows  block  types 
to  be  changed  during  a  run,  but  this  should  not  be  done  unless  there  is  a  good 
physical  reason  for  it;  the  program  may  have  to  be  modified  if  stresses,  for 
example,  must  be  preserved. 

Fully-deformable  blocks  are  subdivided  internally  into  a  mesh  of 
triangular  zones.  Each  zone  is  specified  manually,  using  a  GENERATE  command, 
which  takes  as  arguments  lists  of  grid-points  and  zones.  A  grid-point  must  be 
assigned  to  every  boundary  point  (comer)  of  the  block.  Grid-points  that  are 
placed  on  edges  cause  new  comers  to  be  created. 

4.3  Material  and  Constitutive  Numbers 

Each  block  and  contact  carries  both  a  material  number  and  a  constitutive 
number.  Each  material  number  may  be  associated  with  a  different  set  of  material 
properties.  Independently,  the  program  may  refer  to  different  constitutive  models 
for  different  blocks  or  joints.  The  present  version  of  the  program  contains  block 
and  contact  subroutines  corresponding  to  constitutive  number  1.  These  subroutines 
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model  linear  elasticity  and  Coulomb  friction,  respectively.  Dummy  subroutines 
exist  for  constitutive  numbers  2  to  5;  these  may  be  replaced  by  user-written 
routines  as  explained  in  Appendix  IV.  If  a  contact  is  not  given  a  material 
and  constitutive  number,  it  assumes,  by  default  the  numbers  of  one  of  the  blocks 
comprising  the  contact.  Material  properties  and  even  constitutive  numbers  may 
be  changed  during  a  run,  but  the  user  should  be  sure  that  the  change  is  physically 
reasonable .  The  program  prevents  the  user  changing  block  masses  after  cycling 
has  commenced,  as  this  would  almost  certainly  be  unreasonable. 

4.4  Corner  Rounding 

The  ROUND  command  influences  the  amount  by  which  comers  are  rounded, 
by  setting  the  distance  from  the  true  block  corner  to  the  point  at  which  the 
circle  is  tangent  to  either  edge.  The  command  should  not  be  used  after  cycling 
has  started  because  a  change  in  corner  geometry  would  have  am  unpredictable  effect 
on  contact  forces.  The  rounding  distance  is  also  used  for  other  purposes,  notably 
for  controlling  contact  detection  and  deletion.  Tolerances  are  set  for  these 
functions,  and  in  the  present  version  of  UDEC  are  taken  as  fractions  of  the 
rounding  length.  This  calculation  may  be  changed  by  modifying  subroutine  TOLSET; 
such  a  modification  may  be  necessary  if  very  small  rounding  lengths  are  needed, 
since  in  that  case  contact  updating  may  be  unacceptably  frequent. 

4.5  Printout 

All  input  lines  are  echoed,  and  preceded  by  the  symbol  >  to  distinguish 
them  from  output  produced  by  the  program. 

Printout  is  generally  self-explanatory,  with  full  error  messages. 

f 

4.6  Restrictions  and  Cautions 

The  program  cannot  be  regarded  as  being  in  its  final  form,  as  several 
planned  facilities,  traps  and  options  were  not  completed  in  the  time  available . 

The  program  is  potentially  very  powerful,  as  it  can  model  anything  from  a  con¬ 
tinuum  to  a  complete  discontinuum.  But  this  very  generality  ensures  that  there 
are  many  opportunities  for  misuse . 

The  calculated  critical  time-step  is  only  approximate.  If  numerical 
instability  is  suspected,  the  same  run  should  be  made  with  half  the  time-step 


and  double  the  number  of  cycles.  Any  significant  difference  in  the  results  in¬ 
dicates  that  the  original  time-step  was  too  large. 

Blocks  that  become  detached  from  their  neighbours  may,  under  some  cir¬ 
cumstances,  not  make  correct  contact  again.  The  program  always  maintains  one 
"virtual"  contact  for  a  block  that  becomes  detached;  this  is  to  keep  the  block 
linked  to  the  data  structure.  However  a  potentially  new  contact  that  would 
cross  the  track  of  the  virtual  contact  will  not  be  made  correctly,  if  at  all. 

The  logic  to  deal  with  this  situation  is  quite  straightforward,  but  has  not 
been  written  yet.  No  problems  should  occur  for  fairly  tight  rock  masses. 

At  present,  comer  radii  are  only  re-calculated  when  a  CYCLE  command 
is  given,  although^ the  centre  coordinates  for  each  comer  circle  are  updated  at 
each  time-step.  When  running  problems  in  which  deformable  blocks  are  changing 
shape  rapidly,  the  total  number  of  required  cycles  should  be  split  up  by  using 
several  CYCLE  commands.  In  this  way  the  radii  will  never  be  too  much  in  error. 

4.7  Incompressible  Plastic  Flow 

Fully-deformable  blocks  are  discretised  internally  into  a  mesh  of 
finite-difference  triangles.  Such  assemblies  of  constant-strain  triangles  are 
found  to  be  too  stiff  when  plastic  flow  is  occurring  under  conditions  of  near¬ 
incompressibility:  for  example,  collapse  loads  are  overestimated.  Nagtegaal 
et  al  (1974)  explained  this  phenomenon,  and  Marti  and  Cundall  (1980)  proposed 
a  procedure  for  overcoming  the  problem.  This  procedure  is  called  "mixed  dis¬ 
cretisation"  and  consists  in  averaging  the  volumetric  calculation  over  two 
adjacent  triangles,  while  the  deviatoric  calculation  is  done  separately  for  each 
triangle.  Mixed  discretisation  could  be  incorporated  easily  into  UDEC,  since 
a  linked  list  already  exists  that  could  serve  to  combine  alternate  triangles 
volume trically . 
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5.0  CONCLUSIONS  AND  FUTURE  DEVELOPMENTS 

A  powerful  computer  program  has  been  developed  that  can  model  a  wide 
spectrum  of  problems  in  solid-body  mechanics,  ranging  from  a  continuum  at  one 
extreme  to  a  completely  discontinuous  medium  at  the  other.  Furthermore,  arbi¬ 
trary  mixtures  of  the  two  can  be  accommodated.  The  program  differs  from  pre¬ 
vious  programs  mainly  in  respect  of  its  data  structure,  which  is  designed  to 
have  the  same  topological  properties  as  the  physical  structure  that  it  represents 
The  advantage  in  this  is  that  any  structure  or  connection  that  exists  physically 
has  an  analogue  in  the  linked-list  space  of  the  data  structure. 

The  contractural  objective  of  the  work  reported  here  was  to  develop 
a  program  in  which  rigid  blocks,  simply-deformably  blocks  and  fully-deformable 
blocks  could  interact  with  one  another.  This  was  achieved.  Innovations  that 
have  been  made,  in  addition  to  the  new  data  structure  mentioned  above,  are:  the 
creation  of  block  systems  by  splitting  and  deletion,  and  the  use  of  rounded 
comers  to  prevent  "  locking-up" .  Parametric  studies  may  be  made  in  which  the 

rounding  dimension  is  varied  for  the  same  problem.  More  utilitarian  improvements 
include  free- format  input  with  powerful  parameter  handling  and  continuation  line 
logic,  and  the  fact  that  modifications  and  changes  to  blocks  and  properties  may 
be  made  at  any  stage  during  a  run. 

More  time  was  spent  developing  and  coding  the  new  data  structure  them 
was  anticipated.  Consequently  several  planned  facilities  have  been  provided  for 
but  only  exist  in  skeleton  form  at  present.  The  reason  for  this  is  simply  that 
time  ran  out,  and  not  that  any  difficulty  is  involved.  The  program  UDEC  con¬ 
tains  incomplete  coding  for:  pore-pressure  generation  and  dissipation,  edge/edge 
contact,  automatic  zoning  for  fully-deformable  blocks,  dynamic  cracking,  struc¬ 
tural  connection  and  handling  of  initially-free  blocks .  Some  of  the  coding 
is  almost  complete,  but  in  other  cases  has  only  just  been  started.  The  total 
time  for  completion  of  all  items  noted  above  would  be  about  six  weeks. 

Program  UDEC  has  been  written  in  modular  form  in  anticipation  of  future 
extension.  Further  developments,  such  as  dynamic  input  and  non-reflecting 
boundaries  should  present  little  difficulty. 
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Appendix  1-1 


Input  commands  for  UOEC 


Motes:  Upper-case  letters  in  a  command  or  parameter 

must  be  typed;  the  remaining  letters  are  optional. 
Lower-case  parameters  stand  for  numeric  values.  Integers 
must  be  given  for  parameters  starting  with  l,j,k, l,m,n. 
Real  numbers  may  be  given  as  integers,  but 
not  vi ce  versa. 

Input  is  free-format:  parameters  may  be 
separated  by  any  number  of  the  following 
characters,  in  addition  to  spaces: 

-  C  )  ,  / 

An  additional  line  should  be  given 
at  the  end  of  the  input  file  (after 
the  STOP  command). 

The  first  command  should  be  START  or  RESTART. 


*  *  comment  line 

♦  *  continuation  line 


Block  Material  n  Constitutive  m  xl  yl  x2  y2  ... 

Create  a  rigid  block  of  material  number  n 
and  constitutive  number  m. 

Defaults  are  n*I,  m*  I ,  if  m,  n  omitted. 

Corner  coordinates  are: 

C  x 1 , y 1 ) ,  (x2,y2)  etc., in  a  clockwise  direction.  Continuation 
lines  may  be  used  but  a  pair  of  numbers  defining  a  corner 
must  not  be  separated.  Only  one  BLOCK  command  may 
be  used  per  run  at  present.  Further  blocks  may  be 
created  wi th  the  SPLIT  command,  and  unwanted  ones 
deleted  with  the  DELETE  command. 

Any  blocks  may  be  changed  to  simply-  or  fully- 
deformable  with  the  CHANGE  command. 

Change  xl  x2  yl  y2  Sdef  Material  n  Constitutive  m 

Fdef 

All  blocks  with  centroids  lying  within  the  range 
xl<x<x2  ,  y1<y<y2  are  changed  to  simply-deformable 
or  fully  deformable  (Sdef  or  Fdef  respectively).  Material 
and  constitutive  numbers  may  also  be  changed. 

DAmping  fcrit  freq  Mass 

Stl ff ness 
Internal 

Viscous  damping  is  applied,  in  the  Torm 
of  Rayleigh  damping.  If  a  qualifier  is 
not  given  as  the  third  parameter,  full 
damping  is  used,  with  fcrit  as  the  fraction 
of  critical  damping,  and  freq  as  the  centre 
frequency.  The  word  "Mass"  eliminates  the 
stiffness-proportional  dashpots,  and  "Stiffness" 
eliminates  the  mass-proportional  dashpots. 
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The  word  "Internal"  causes  the  specified 
damping  to  be  applied  to  the  3  internal  degrees 
of  freedom  of  s impl y-deformabl e  blocks. 

DEIete  xl  x2  yl  y2 

All  blocks  are  deleted  in  the  range 
x1<x<x2,  y 1 <y<y2 

Dump  n  m 

Dump  memory  to  printer  from  the  main  array 
from  address  n  to  addess  m.  Internal 
pointers  MFREE,  JUNK,  IBPNT,  ICPNT  and 
IDPNT  are  also  printed.  MFREE  gives  the 
highest  memory  location  that  is  currently 
free. 

FRAction  f 

f  is  taken  as  the  fraction  of  critical 
time-step  to  be  used. 

Fix  xl  x2  yl  y2 

All  blocks  are  fixed  in  range 
x1<x<x2  ,  yl<y<y2 

FRee  xl  x2  yl  y2 

All  blocks  are  set  free  in  range 
xl<x<x2  ,  yl<y<y2. 

Note:  by  default,  all  blocks  are  free  initially. 

Generate  xl  x2  yl  y2  Manual  Gridpoints  <glist>  Zones  <zlist> 

Au  toma  t i c 

The  first  block  encountered  in  the  range  xl<x<x2,  y!<y<y2 
is  discretised  as  ful ly-deformabl e.  The  automatic 
option  is  not  available  yet.  For  manual  generation, 
a  list  of  grid-points,  <glist>,  and  zones,  <zlist> 
must  be  given.  The  format  for  <gl!st>  is: 

xl  yl  x2  y2  x3  y3  ....  , 
where  each  x,y  pair  is  a  coordinate  of  a  grid-point. 

If  a  given  coordinate  lies  within  a  certain 
tolerance  of  a  block  corner,  the  grid-point  is  placed 
on  that  corner.  If  the  coordinate  lies  within  the 
same  tolerance  of  a  block  edge,  a  new  corner  is  created  In 
the  edge.  The  tolerance  is  taken  as  0.9  times  the 
rounding  length.  The  format  for  <zllst>  is: 

1 1  ml  nl  1 2  m2  n2  . 

Each  triple  corresponds  to  the  three  grid-points  that 
define  the  zone,  where  the  numbering  of  the  grid- 
points  refers  to  the  order  in  <g1ist>,  starting  with 
the  last  (l.e.  last  grid-point  Is  number  1). 

Both  <glist>  and  <zllst>  may  extend  over  an 
arbi trary  nunber  of  continuation  lines,  but  doubles 
and  triples  should  not  be  spilt  over  two  lines. 
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Gravi ty  gx  gy 

Gravi tational  accelerations  are  set  for 
the  x-  and  y-  directions. 

°Lot 

All  blocks  and  centroids  are  plotted 

Print  Blocks  Contacts  CORners  Domains  List  DLlst 

Data  are  printed  on  blocks,  contacts,  corners, 
domains  and  linked  lists  for  blocks  and  domains. 


PROperty  Material  n 
n 


Sulk  b  Cohesion  c 
K  »b  Dens  I ty  d 
KN»sn  Friction  f 
KS«ss  G*g 


Material  properties  are  defined  for 
material  number  n.  Properties  are: 
bulk  modulus, b;  shear  modulus, g;  density, d; 
joint  normal  stiff ness, sn;  shear  st i f f ness, ss; 
friction  coeff i ci ent, f ;  cohesion, c. 

The  first  parameter  must  be  the  specification 
of  material  number. 


Restart 

The  program  is  restarted,  using  data 
from  the  restart  file 

RSet  v  ia  ioff 


The  real  value  v  is  inserted  in  the 

main  array  at  address  fa,  with  offset  ioff. 

ROund  d 

Each  block  corner  is  rounded  with  a 
circle  that  is  tangential  to  the  two 
correspond! ng  edges  at  a  distance  d 
from  the  corner. 

SAve 

The  current  problem  state  is  saved  on 
the  restart  file. 

SCale  s 


Plot  scale  is  set  to  s 

SPlit  x 1  yl  x2  y2 

All  blocks  in  the  path  of  a  line 
extending  from  point  (xi  yl )  to  (x2  y2) 
are  spilt  Into  two.  At  present,  the  line 
should  not  pass  through  any  corner,  or 
run  too  close  to  an  existing  edge. 
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STArt 


The  program  does  a  cold  start. 
Stop 

The  run  stops. 
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Appendix  1-2 


Error  numbers 


1  Memory  overflow 

2  Unrecognisable  command 

3  Mot  start  or  restart  as  first  command 

4  Material  number  out  of  range 

5  Block  has  less  than  3  corners 

6  Negative  or  zero  block  area 

7  Missing  parameter 

8  Missi ng  y-value 

9  Constitutive  number  out  of  range 

10  Unrecognised  parameter 

11  Contact  stiffnesses  undefined  -  cannot  cycle 

12  Zero  mass  block(s)  present  -  cannot  cycle 

13  Mass  damping  for  rigid-body  motion  too  high. 

14  Rounding  length  too  great 

15  Contact  overlap  too  great. 

16  Internal  mass  damping  too  high. 

17  Cannot  split  fu 11y-def ormab 1 e  block. 

18  Only  one  block  may  be  created  at  present  - 

use  SPLI T  for  more. 

19  Cannot  delete  final  contact  in  problem. 

20  Internal  error,  subroutine  DELC. 

21  Cannot  delete  ful ly-deformable  block  at  present. 

22  Fifth  parameter  must  be  "Manual"  or  "Automatic" 

23  Not  available  yet 

24  Cannot  find  block  in  range 

25  Zone  pointer  references  non-existent  grid-point 

26  Mi ss i ng  data 
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APPENDIX  II  EXAMPLE  RUNS 


] 

j 


II-l  Ball  Rolling 

II-2  Collapse  of  Opening  in  Jointed  Rock 


This  example  illustrates  that  corner  rounding  can  be  taken  to  the 
extreme  of  creating  a  circle  from  an  angular  block.  The  complete  input  se¬ 
quence  to  create  the  blocks  and  run  the  problem  is  given  overleaf.  Fig.  II- 1 
shows  the  system  of  blocks  after  splitting,  but  before  deletion  of  unwanted 
blocks  and  final  rounding.  Fig.  II-2  is  the  plotted  output  from  UEEC,  with 
plots  superimposed  after  successive  200-cycle  increments. 


S  f  A  rt  T 

PHOP  MAT  =  1  DENS=2000  KN=i.C08  KS=t.E08  F=2. 

A V  I T V  *.  -8. 

DAMP  =.5  16.  MASS 

PH  AC  0.10 

uLuCK  0. i *2.  0  •  ,  1 0  •  5  14., 10. 5  1 4 « / *2. 

SPL 1 l  -1. ,  .5  lS.  ,8.5 
SPLIT  10.5, S.  7.5,11. 

SPLIT  12.5,6.  10. ,11. 

SPLIT  -l . ,3.  12.  ,9.5 

Dt.Le.re,  0.,9.  4>.  ,  7  . 

bfc.Lt l f.  0.  ,9.  6.  ,  10. 

DlLEIc  8.  ,  H  .  8.  ,  1  !  . 

OtLETt  U.,14.  7. ,11. 

e  I  a  u .  ,  1 4 .  0.  ,  6. 

kOUijD  1.115 
C'ICLt  300 
P  n 

PLOT  0 

PKUP  MAisl  r=0.8o 
DAMP  0.1  3>,‘  S  i  Lr  F'lh^ 

GKAVili  0.  -10. 

CTCbt  200 
P  D 

PLOT  1 
CTCLt  2oo 

P  D 

PLOT  1 
CTCbt  2'>ij 
P  D 

PLOT  1 
CiCLc  2uu 
P  o 
PLOT  1 
CYCLE  tuO 
P  6 

plut  i 

SAVE 

STOP 

END 
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FIG.  II-I!  SYSTEM  OF  BLOCKS  PRIOR  TO  DELETION 


FIG. II -2=  MOTION  OF  'BLOCK' UNOER  GRAVITY 


II-2  Collapse  of  Opening  in  Jointed  Rock 

Figure  II-3  shows  the  initial  state  of  a  rock  system  as  produced  by 
the  input  sequence  given  on  the  following  page.  The  example  is  intended  to 
demonstrate  the  interaction  of  all  three  block  types  in  the  same  problem.  T 
upper  layer  of  rigid  blocks  is  heavier  than  the  rest.  Figures  II-4  onwards 
record  the  movements  that  take  place  after  several  blocks  have  been  removed. 
It  is  interesting  to  note  that  "hang-ups"  do  not  occur,  due  to  the  rounded 
corners.  The  fully-deformable  blocks,  which  are  discretised  internally  into 


START 

SCALE 

*  SMALu  tLt>l  Uf  RilvAij 
ROUND  1.3 


FRAC  .030 
GRAVlii  j  -It 

BLOCH  o  0  '<  40  cO  -t'J  cu  t 

PROP  i-i a'I  =  L  i-i.  ijitO'jO  iViSl.Or.h  i\b5  1.i<t.e  t;U iiK  —  1  •  utV  U  —  0.4e.l  F 
PROP  !«iasi  UiV>S=5tuO  Kt.sl  .Ut0  fxa=l  .UC.I;  h=U.3 


SPLIT  -) 

10 

o  1  It 

SPLIT  -l 

z  0 

o  1  2  0 

SPLIT  -1 

30 

3 1  3') 

SPLIT  lu 

4 

It  Hi 

SPLIT  20 

■3 

20  4 1 

SPLIT  it 

'l 

30  4 1 

SPLIT  40 

'I 

4:J  4  1 

SPLIT  bo 

V 

30  i  1 

SPLIT  20 

7.0 

60  4  3 

.  b 

CHANGE  0 

CM* 

30  40 

m A  TER t  k  t—  2 

CHANGE  o 

3  0 

20  30 

SDEr 

Change  o 

20 

It  2o 

FDEF 

CHANGE  2o  *0  IJ  z u  Suc.r 
CHANGE  4U  ou  lu  2‘i  FOcF 

FIX  0  o  0  o  It 


GENERATE  t  it  lu  zt  At  GKiOPuXwiS  (.3*15)  (t,/.0)  (0,10) 

♦  (It, It)  lit, zt)  ZuoES  1  *  x ,  b  2,3,3  3*4,5  *,1,3 

GENERAl't.  1  o  '4  0  It  40  ••  G  lb,  lb  lu,zt  lu,lu  20,  It 

^  20  *  2t  TJi^imiS  1*2*5  2*3*3  3*4*3  4  *  1  *  b 


GENERATE  <*0  30  10  20  rt  t  4b,  lb  4t*/u  40,10 

^  bO  *  z  1*4*5  4*3*3  3  *  4  *  S  4*1*3 

GENERATE  oO.oO  10,20  .*  o  35*13  5t*zi.  30,10 
4  50*20  Z  1*2*5  2*3*5  3*4*5  4*1*3 


DUMP  1  1 

DAMP  .25  , 
DAMP  0.25 
CYCLE  l'>o 
PLOT  0 
P  B 
STOP 
END 


1.0  1  i  i  L  R !'  A  L 


7*^ 


i  V  Xi ' 

i 


30,10 
o  0 , 1  u 
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APPENDIX  III  -  PROGRAM  GUIDE 


III-l 

III-2 

III-3 

III-4 


Main  Common  Block  Variables 
Parameters  &  Data  Groups 
Subroutine  Functions 
Subroutine  Calling  Map 


50  - 


LI  NEC  80)  Buffer  for  current  input  line  in  A1  format. 
LINEK80)  Buffer  for  next  input  line. 

LPNT( I )  Pointer  to  start  of  parameter  I  in  L!NE(  ) 
after  removal  of  blanks,  etc. 

RAFLAG 

PPF LAG  .TRUE,  if  pore-pressure  calculation  requested 

ERFLAG  .TRUE,  if  an  error  has  occured 

STFLAG  .TRUE,  if  the  first  input  line  has  been  processed 

COFLAG  .TRUE,  if  the  current  line  is  a  continuation 

NCFLAG  .TRUE,  if  the  next  line  is  a  continuation 

JMPSAV  Index  of  last  computed  GOTO  in  MOM 

NERR  Error  number 

JUNK  Pointer  to  list  of  spare  memory  groups 

MFREE  First  unused  memory  address 

I  BLOCK  Current  block  number 

I  DOM  Current  domain  number. 

I  STACK  Stack  pointer 

NCYC  Currently  requested  number  of  cycles 

MCTOT  Total  number  of  cycles 

TDEL  Time-step 

FRAC  Requested  fraction  of  critical  time-step 

I  ROUTE  Routing  number,  used  in  UDC 

NLINE  Output  line  count 

NPAGE  Output  page  count 

JMPGEN  Routing  number  for  continuation  line  in  GEN 
ALPHA  Mass  damping  coefficient 

BETA  Stiffness  damping  coeffient 

CONI  Damping  factor 

C0N2  Damping  factor 

ALPB  Internal  mass  damping  coefficient  for 

simply  deformable  blocks. 

C1B  Damping  factor  derived  from  ALPB. 

C2B  "  " 

DEGRAD  PI/ 180 

PI  3.14159 

PSCALE  Plotting  scale 

ATOL  Distance  between  particles  at  which  a  contact 

is  first  formed. 

BTOL  Distance  between  particles  at  which  a  contact 

is  broken, 

CTOL  Maximum  (negative)  overlap  allowed 

when  forming  contacts 
DTOL  Rounding  length 

DT0L2  «0TOL/2.0  (maximum  contact  overlap) 

ETOL  Limit  on  maximum  domain  displacement 

to  trigger  contact  update. 

FTOL 

GTOL 

HTOL 
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I BPMT  Pointer  to  list  of  blocks 

ICPNT  Pointer  to  list  of  contacts 

IDPNT  Pointer  to  list  of  domains 

I00PMT  Pointer  to  outer  domain 

AKN( I )  Normal  joint  stiffness#  material  I 

AKS(  I  )  Shear  joint  stiffness#  material  I 

AMU( I  )  Joint  friction  coefficient#  material  1 

C0H( I )  Joint  cohesion#  material  1 

0ENS( I )  Density#  material  I 

BULK(I)  Bulk  modulus#  material  I 

SHEAR(I)  Shear  modulus#  material  I 

A  LAM  1(1)  Lame  constant#  material  I 

ALAM2 ( I )  Lame  constant,  material  i 

l A(  )  Mai n  array 


-  52 


f 


Appendix  III-2  Parameters  a  Data  Group 


Offsets  for  block  data  array 


Note:  the  first  integer  in  each  block  array 
— —  (offset  0)  is  the  block  type  number,  as  follows 

1  rigid  block 

2  s imp! y-deformabl e  block 

3  fu 1 1 y-deformabl e  block 

KB  Pointer  to  next  block  in  block  list 

KP  Pointer  to  one  corner  in  block's  corner  list 

KMAT  Material  number 

KCOMS  Constitutive  number 

KX  x  coordinate  of  centroid 

KY  y  coordinate  of  centroid 

KXP  x  veloc i ty 

KYD  y  veloci ty 

KTD  Angular  velocity  (anticlockwise  positive) 

KSM  Block  mass 

KB  I  Moment  of  inertia 

KBFX  x  centroid  force-sum 

K8FY  y  centroid  force-sum 

KBFT  centroid  moment  sum 

K8EX  extension  pointer  (to  SDEF  or  FDEF  data) 

Offsets  for  corner  data  array 


Note:  the  first  integer  (offset  0)  contains 
----  the  value  MCOR  to  denote  a  corner 

KL  Pointer  to  next  corner  or  contact  on 
block,  in  clockwise  direction. 

KR  Pointer  to  next  corner  in  anticlockwise 
dl rection 

KNB  Pointer  to  host  block 

KXP  x  coordinate  of  corner 

KYP  y  coordinate  of  corner 

KXCP  x  coordinate  of  local  circle  centre 

KYCP  y  coordinate  of  local  circle  centre 

KRAD  Radius  of  local  circle 

KXDP  x  veloci ty  of  corner 

KYOP  y  veloci ty  of  corner 

KGP  Pointer  to  corresponding  grid-point  if  block 
is  fu 1 ly-deformab le 
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Offsets  for  contact  data  array 


Mote:  the  first  integer  (offset  0>  contains 
the  value  MCON,  to  denote  a  contact 

KC  Pointer  to  next  contact  in  contact  list 

K31  Address  of  first  block  Involved  in  contact 

KB2  Address  of  second  block  Involved  in  contact 

KL1  Pointer  to  next  1  ten  in  clockwise  list 

of  block  correspond! ng  to  KB  1 
KL2  sane  as  KL1,  but  for  block  KB2 

KD 1  Address  of  domain  to  left  of  contact, 

going  from  block  KB1  to  KB2 
KD2  Address  of  domain  to  right  of  contact 
going  from  block  KB1  to  KB2 
XCM  Material  type  number 

KCC  Constitutive  number 

KXC  x  contact  coordinate 

KYC  y  contact  coordinate 

KXDC  relative  x  velocity  (of  block  KB2  relative 
to  block  K81) 

KYDC  relative  y  velocity 

KCS  relative  shear  displacement 

KCN  relative  normal  displacement 

KCFS  shear  force 

KCFN  normal  force  (positive  compression) 

XCCOD  code  number: 

1  corner /corner  contact 

2  corner/edge  contact  ( KB  1  ..  corner, 

KB2  ..  edge) 

3  edge/corner  contact 
Offsets  for  domain  data  array 


Note:  the  first  integer  (offset  0)  contains  the 
----  value  MDOM,  to  denote  a  domain 

KD  Pointer  to  next  domain  in  domain  list 
KPP  Pore-pressure  for  domain 
KUMAX  Fictitious  domain  displacement 
KDL00P  Pointer  to  one  contact  in  anticlockwise 
list  around  domain 


-  54  - 


Simply-deformable  extension  array 


K  ED  1  1 

) 

KED12 

) 

S  trai n-rate 

KED21 

) 

tensor 

XED22 

) 

ksiii 

) 

KSI  12 

) 

Internal  stress 

KSI  21 

) 

tensor 

KSI  22 

) 

KSA1 1 

) 

Applied  stress 

KSA12 

) 

tensor  (multiplied 

K3A21 

) 

by  block 

KSA22 

) 

area) 

Offsets  for  grid-point  data 


KG  Pointer  to  next  grid-point  in  grid-point  list 

KCOR  Pointer  to  corresponding  block  corner 

KXG  x  coordinate 

KYG  y  coordinate 

KXDG  x  veloci ty 

KYOG  y  veloci ty 

KGFX  x  force-sum 

KGFY  y  force-sum 

KG PM  grid-point  mass 

Offsets  for  zone  data 


K2  Pointer  to  next  zone  in  zone  list 

KZG  Start  of  triple  pointer  to  3  surrounding 

gr i d-ooi nts 
KZS11  ) 

KZS12  )  Stress  tensor 
KZS22  ) 

KZM  Zone  mass 

Logical  unit  numbers 


LUNIF  Unit  number  for  input  file 

LUNOF  Unit  number  for  output  file 

LUNG  Unit  number  for  general  I/O  (e.g.  restart) 

LUMP  Unit  number  for  plotted  output 


Number  of  words  In  data  arrays 


NVCR 

Corner 

MVBL 

Block 

MVCN 

Contact 

MVDO 

Doma i n 

extension 

NVSD 

Simplydeformab  le 

NVZO 

Zone 

NVGP 

Gr  id-poi nt 

Array  limits 

MTOP 

Size  of  main  array  ( ! A) 

NMAT 

Maximum  number  of 

mater i als 

NCONS  Maximum  constitutive  numbers 

NTYP  Number  of  block  types  (rigid,  SDEF  etc) 

Head  codes  (contents  of  first  integer  in  data  array) 


•  m 

MR  1  G 

*1  Rigid  block 

MSDEF 

*2  Simply-deformable  block 

mfdef 

*3  Fu 1 ly-deformab le  block 

MCOR 

Corner 

MCON 

Contact 

MOOM 

Doma i n 

l 


sndix  III-3  Subroutine  Functions 


SUBROUTINE  ACCOM ( A ) 

-RECuKD  MAX  AND  M I M  VELOCITIES  fuK  CORNER 


SUBROUTINE  aPLGT 

-BLOCKS  AkE  PLOTTED 


SUBROUTINE  akCIRAOJ 

-PLOT  ARC,  kADIUo  KAO ,  CENTkE  (Ap,itO,  FROM  (XI, Yi) 
-TO  (X2,Y2J.  PEN  ASSUMED  du*  N  ON  e.NTKi. 


SLOCK  DATA  FrtfcD 

-INITIALISE  FIXED  PARAMETERS,  iHlo  ROUTINE 
-(AND  THE  COMMON  SLOCK  /PAkAM/.J  CAn  6c.  REPLACED 
-BY  A  SET  ur  PARAMETER  STATEMENTS  IF 
-SUPPORTED  a Y  THE  COMPILER. 


SUBROUTINE  BRIDGE  ( IAD 1  , 1  AD  2, 1  AST  ,  Iao2 , ICOL , IDNE* ) 

-MAKE  A  CONTACT  TO  BRIDLE  DOMAIN  IUuM  BETmEEN  LIST  ITEMS 
-IA01  AND  I AD’2  (BLOCKS  IASI  M  iAti2J.  ICUU=C0D£  OF  N£w  SEC, 
- (XC , YCJ sCOQRDlNATtS  OF  CONTACT. 

-IDNENSRETURNED  n£I«  DOMAIN  NUrtoEn. 


SUBROUTINE  CCCA1,.A2) 

-PARAMETERS  FOK  CORNER-CORNER  CONTACT 


SUBROUTINE  CECA1,A2,AP) 

-PARAMETERS  FOR  CORnER-EDGE  CONTACT 


SUBROUTINE  ChANGE 

-CHANGe.  ATTRIBUTES  OF  EXISTING  BLuCrvS  allnlm  RANGE 


SUBROUTINE  CL1 

•ELASTIC,  ISOTROPIC  CONSTITUTIVE  LA« 


i 

ss 
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SUBROUTINE  CQORO(IAC) 

-UPDATE  COukLInmTeS  FOR  CONTACT  iAC 


SUBROUTINE  CKAD(IAB) 

•CORNER  RADIUS  AND  CIRCLE  CENihfc  row  EACn 
-CORNER  UF  SLOCK  IAo 


SUBROUTINE  CKfcAlt 

-CREATE  A  o l, U C A  (FROM  dLQCk  CummAnuj 


SUBROUTINE  CVEL(B,C) 

-UPDATE  RELATIVE  CONTACT  VELOCITIES 

-B(  )  IS  BLOCK  ARRAY;  C(  )  IS  CONTACT  ARRAY 
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subroutine;  cvelfdccj 

-RELATIVE  CONTACT  '/ELuCITIe.S  From  G.P.  VELOCITIES 


SUBROUTINE  CYCLe. 

-MAIN  CALCULATION  CYCLE 


SUBROUTINE.  DAi-lP 

-PKOCfcSS  DAMPING  COMMAND 


SUBROUTINE  DElB 

-DELETE  BLOCK  IBLOCK 


SUBROUTINE  DELCtlAC) 
-DELETE  CONTACT  IAC 


jtg  ' 

vX  pG 


SUBROUTINE  DELLSTCIAD, IPNT,KQFF) 

-DELETE  ITEM  IAD  FROM  LIST  POINTED  TO 

-BY  IPNT  ,  WITH  OFFSET  WITHIN  AN  Il&M  OF  KuFF 


SUBROUTINE  FDCCIPP) 

-DEAL  nITH  FORCE  APPLIED  AT  CuRNtK  OF  FULLY -DEF  BLOCr 


SUBROUTINE  FOECIPA,IPB) 

-DEAL  WITH  FORCE  APPLIED  TO  EDGE  OF  FULL Y-DEF  BlOCK 


SUBROUTINE  FDEF 

-SCAN  THROUGH  GRID-POINTS  &  ZUNtS  FUK  FULLY-OEF  BLOCK 


SUBROUTINE  FINDCN  ,  NG) 

-FIND  N  WORDS  OF  MEMORY 
-NG  -  RETURNED  ADDRESS 

-ERFLAG  IS  SET  IF  MEMORY  CANNOT  BE  FOUND 


SUBROUTINE  FIX(IFIX) 

-SET  OR  RESET  FIX  FLAG  FUR  ALL  bLUCKS  WITHIN  RANGE 


SUBROUTINE  FORDCCON) 

-FORCE/DISP  CALC.  FOR  CONTACT 


SUBROUTINE  GEN 

-PROCESS  "GENERATE"  COMMAND,  wnlCH  GENERATES  ZONING  FOR 
-FULLY-OEFOkMABLE  BLOCK  EITHER  aAnuaLLY  Ok  AUTOMATICALLY. 


SUBROUTINE  Gh.0  C I AB  ,  Ak  ,  KG2  ) 

-FUK  BLuCK  and,  COMPUTES  AREA,  RAuYUo  OF 
-GYRATION  SUUARED  ANO  CENTROID  tAP.iPJ 


FUNCTION  GETR(A,KOFF) 

-GET  REAL  VALUE  FROM  MAIN  ARRAY  wITH 
-OFFSET  KOFF 
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subroutine:  halt 

-EKRGk  MESSAGES  PRINTED  HERE 

-PROGRAM  STOPS  UNDER  CERTAIN  CuNDlTIUi'.S 


FUNCTION  ICALo(IC) 

-GET  CALLING  ADDRESS  FOR  CONTACT  IC  FROM  SLUCn  IBLUCK 


FUNCTION  IDUmP(IAD) 

-RETuKN  DOMAIN  ADDRESS  FOR  SEGMENT  I AD  ,  cLOCn  IbLLCR 


FUNCTION  IGE’I  C IAD , KuFF ) 

-GET  INTEGER  FROM  MAIN  ARRAY  A  lAu+EUrr 


FUNCTION  1GETGPCI) 

-FIND  I'TH  GRIDPOINT  ADDRESS,  BLUCn  IBLOCK 


SUBROUTINE  IMSG 

-INITIAL  MESSAGE  FROM  CODE 


SUBROUTINE  INERT C XU , YO , X , Y , YP , YS , kmIL’aa ) 

-  COMPUTE  MOMENT  OF  INERTIA  AouUT  GIVEN  AXIS 


SUBROUTINE  INI 

-PREPARE  FOR  CYCLING 


subroutine  insect ( yes ) 

-find  IF  2  LINES  INTERSECT  CkETURwS  YES*. TRUE.) 
-AND  IF  SO,  RETURN  COORDS  (XP,xP). 

-THE  2  LINES  ARE  C (XI , Y2 ) , CX2 , Y2 J ) 

-  AND  ( CX3 ,Y3),(A4,Y4)) 


SUBROUTINE  10(1, NPAR) 

-ALL  NON-STANDARD  I/O  DONE  HERE 


FUNCTION  IVARCNPAR) 

-TO  RETURN  INTEGER  VALUE  OF  PARAMETER  nPAR 


SUBROUTINE  JUMPS 

-TO  JUMP  TO  NE#  BLuCK  IN  DOMAIN  FRuM  CONTACT  I AD 


SUBROUTINE  LOS£(N,NG) 

-TO  RETURN  N  ROKOS  OF  MEMORY  AT  AdDKESS  NG 


SUBROUTINE  MAREo(IAb,UKAS,UM01,A, X ) 
•Ttl  Tm&IaLL  MASS,  mijI  o  CewIkuId 
-IN  BLOCK  I AB 


SUBROUTINE  MAKECC IAC  ,IAfii_/IAb2J  lul^X-Vix  IDlj  ID  2^,  I  COD) 

•FILL  CONTACT  DATA  BLOCK  AT  IAC,  AND  LINK  TO  CONTACT  LIST. 

— IAB1 ,  IAa2  s  BLOCK  ADDRESSES 

-IL1  ,  IL2  3  ASSOCIATED  CIRCULATING  LISTS 

-ID1  ,  ID2  s  DOMAIN  ADRESSES  -  59  - 

-ICUD  3  CONDITION  COOE 


SUBROUTINE  MAK£P(IAO,IAL,IaR,X, 

>FORH  CUKNeiK  0  I  AD ,  WHERE  1  AL=CLOCkw  ISE  LINK, 
-IARsREVtRSE  LInK,  C X, K)=CORNEh  CDuRuS, 
-IA6=bUUCK  ADDRESS 


SUBROUTINE  MATCH  (  NTAB  ,  INOEX  ,  NPAK  ,  DUMP  ,  oAu  ) 

-TO  MATCn  INPUT  SThI.mG  TO  KEiwukd  It*  TABLa. 
-INPUT:  NTmB  TABLE  OF  KEYWORDS 

-  INDEX  LIST  Of  LtNUTHS 

-  NR AR  PARAMETER  NO.  IN  INRUT  LINE 
-OUTPUT:  JUMP  DISPATCH  NUMBER 

-  BAD  .TRUE.  FUR  MISS I No  PARAMETER 

-  OR  STRING  NOT  FOUND. 


SUBROUTINE  MUN 
-MONITOR 


SUBROUTINE  MOTIQN(A) 

-LAw  OF  MOTION  FOR  SINGLE  BLuCK 


■,T 

'•y  ■ 

£4 


Of  £ 


SUBROUTINE  MOVE ( IAD 1 , I AD2 , NWO 3 

-MOVE  NwD  WORDS  FROM  ADDRESS  IAD1  TO  IAD2 


SUBROUTINE  MVBAK(IAC,IAB,ICQO) 

-MOVE  CONTACT  BACKWARDS  PAST  CORNER  (TAKING 
-WITH  IT  ALL  INTERMEDIATE  COw TACTS). 
-lACsCUNTACT;  IABsbLOCK;  ICOUsnLw  CODE. 


SUBROUTINE  MVFUR(IAC,IA6,IC00) 

-MOVE  CONTACT  FORWARDS  PAST  CuKNtr  (TAKING  wITH  IT 
-ALL  INTERMEDIATE  CONTACTS) . 

-lACsCUNTACT;  I AB-bLOCK ;  ICOu  =  nEw  CUDc. 


SUBROUTINE  NEXTCP (IAD, IADN, ICALL ) 

-FINO  NEXT  ITEM  ( IADN )  IN  CORNER  CriAlN  AND 
•ITS  CALLING  ADDRESS,  GIVEN  CURRENT  ITEM,  IAD 
-BLOCK  NUMBER  ASSUMED  TO  BE  IblUCK 


SUBROUTINE  NEXTDCIAD, IADN) 

-FINO  NEXT  ITEM  (IADN)  IN  DOMAiN  CHAIN. 
-IADsCURRENT  ITEM;  IdOMsDOMAIn. 

-IBLOCK  IS  UPDATED  FOR  EACH  CONTACT  ENCOUNTERED. 


SUBROUTINE.  nEXTuS 

-GET  NEXT  SEGMENT  IN  DOMAIN,  «ntnt  oe.GW'tuT  IS 
-DEFINED  AS  CORNER  OR  EDGE.  ImS  1b  oc.Gnc.nT  NuMscR. 
-CORNERS ,  IrUc. .  IF  CURRENT  SeGmbM'  IS  A  CuKncR, 
-FALSE  OTHERWISE.  CONXAC*.TRUt.  If  A  CONTACT 
-HAS  JUST  BEEN  PASSED  THROUGH. 

-"ITEM"  COUNTS  CONTACTS+SEGMENTS 
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SUBROUTING  Nc.xrP(IAO,IAP) 

-FIND  NeXT  CCKNEk  CIA?)  In  BLOCK  CHAIN, 
-GIVEN  CuKHc.Nl’  ITEm,  IAO,  IN  CHAIN. 
-BLOCk  nuMoEk  IS  IoLOCk 


SUBROUTING  OUTCP(IAO) 
-PRINT  LIST  ITEM  IAO 


SUBROUTING  UUTruGUADJU 

-SCAn  THKGuGn  GRID-PGIMS  FOR  PHIi.luUl 


SUBROUTINE  GUTFOZ ( I AO  1 ) 

-SCAN  THROUGH  ZGnE  FOR  PRINTOUT 


SUBROUTINE  OU TG ( A , I AO , ICQR ) 

-PRINT  DATA  FOR  GRIO-POInT,  akkai  aO 


SUBROUTINE  OUTSO  C  A ) 

-PRINT  OUT  INTERNAL  VARIABLES  t OR  SIMPLY-CLF.  bLGCR 


SUBROUTINE  ?LOTS(X,Y,I) 
-SCALEO  PLOT 


SUBROUTINE  POP (ITEM) 

-POP  ITEM  FROM  STACK 


SUBROUTINE  PRINT 

-RAIN  PRINTOUT  IS  DONE  HERE 


SUBROUTINE  PROP 

-RESPUNU  Tu  •'PROPERTY"  COMMAND 


SUBROUTINE  PUTFCEC8) 

-PUT  FORCES  r  RUM  CONTACT  INTO  oLUCK  AT  Bl  ) 


SUBROUTINE  PUTFDF(A,FX,FY) 

-AUO  FORCES  INTO  GRIO-POInT' S  GOhCt  SUM 


SUBROUTINE  PUTI(IAO,KOFF, I) 

-PUT  INTEGER  I  AT  IACIAO+KOFF) 


SUBROUTINE  Push (ITEM) 
-PUSH  litrt  JNTU  STACK 


SUBROUTING.  Pu  rptAKKAi  ,  I ) 

-PUT  VALUE  OF  PARAMETER  I  In  akKaY 
-ANO  CHECK  FOR  RANGE  OF  I 


SUBRflUT IN E_  PUTRiA^KOFF  t  V ) 

-INSERT  HEAL  VALUE  V  InTO  MAIN  ARRAY 
-WITH  OFFSET  KOFF 


SUBRQuTINt  PuTSTRCS) 

-PUT  APPLIEj  STRESSES  IN  SIMPLY -OcFUKMABLL  aLOCis 


SUBROUTINE  PVt uCo,C) 

-UPDATl  CUK.iEK  VELOCITIES  AND  CUOruIN-TES. 
-BC  )  IS  BLOCK  ARRAY;  CC  )  IS  Cuhnc.H  Akka  i . 


SUBROUTINE  PVELFDCm) 

-UPuATE  CORNER  CuGKOS  &  VELGCii'iEO  tsun  L.f.  VmLUlo 


LOGICAL  FUNCTION  QoCCIAb) 

-RETURNS  .TkUE.  IF  BLOCK  IAB  nAS  ANY  CUMaCTb 


SUBROUTINE  QUERY ( IAC ) 

-CORNER  IS  SET  .TRUE.  IF  CONTACT  IAC  ON 
-BLOCK  ISLOCK  CORRESPONDS  TO  A  CORNER. 


LOGICAL  FUNCTION  SANGECIAB) 

-.TRUE.  IF  nLOCK  CENTROID  LI £.5  OuIbIUE  RANGE  LIVgN  til 
-USER. 


SUBROUTINE  KEL1 

-COORDS  Or  POINT  CAP,YP)  RELATIVE  TO 
-LINE  C CXI, fl) , CX2, Y2) ) .  COMPUTES  ALSO 
-LENGTH  CZ),  COS  CC)  &  SIN  CS)  Of  LlNc. 
-OUTPUT  COORDS;  CAR, YR) 


SUBROUTINE  REL2 

-SAME  AS  RELl,  BUT  C  AnU  S  ARE  ASSUMED  AS  INPUT 


FUNCTION  RVARUPAR) 

-TO  RETurN  REAL  VALUE  OF  PakamEIer  n PAR 


LOGICAL  FUNCTION  SEPCICHAR) 

-RETURNS  .TRUE.  IF  ICHAR  IS  A  SEPARATOR 


SUBROUTINE  SPLIT 

-SPLIT  BLuCkS  In  THE  PATH  OF  GIVEN  LINE 


SUBROUTINE  SPLT C I ADI , IAD2 , XA , YA , Xd , YB ) 

•SPLIT  BLOCK  CIBLOCK)  BETWEEN  Um,  1a)  AND  CAB, Y  B ) 
•IAD1  AND  IAU2  ARE  LIST  ITEMS  uuSL  beFORE 
-SEGMENTS  inAT  .»Ilu  BE  SPLIT 


SUBRUUL'Inc.  oImrI 

•STARTUP  ASSIGNMENTS  i  INITIALISATION 


SUBROUTINE  STKAINCB,e) 

•COMPUTE  STRAIN-RATE  TENSOR  FOR  SIMPL Y-Dtf ORMABlE  bLuCK 
-BC)  IS  THE  BLOCK  ARRAY 

-EC)  IS  The  EXTENSION  ARRAY  Fur  SlMr-LY-DEF  DATA 
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SUBROUTINE  STke.SS(.B,£) 

-STKeSSeS  F  rtuM  STRAIN-RATES  ...  olMPbY-Dci'  (.•RmAeLE 

-BO  IS  THe  BLuC.K  ARRAY 

-£()  IS  Trie  SDEr  EXTENSION  AkRax 


oLOCK 


SUBROUTINE  TIu>Y  . 

-Xu  c.Ll.*i  I  hi*  i’c.  o LANK'S#  ire.  t  Kom  lue’Ur 

-LI.^E  and  InOea  TO  bUCAliuN  UF  R  AK  A  MET  E  h  S 


SUBROUTINE  TOLSET 

-S£T  l'Ulibri AiiCES  FROM  AVEkaGE  auOCK  DIMENSIONS 


SUBROUTINE  TRUE ( X AO , X  ,  i  ) 

-GET  TRUE  blST  ITEM  BEFORE  CUT  a  (X,Y). 
-IAl)  IS  WEAKEST  CURnEk  oErURE  (A,Y),  AnL 
-IAU  IS  OVERWRITTEN  wITh  New  LIST  AOukESS. 


PROGRAM  uDEC 

-UNIVERSAL  DISTINCT  ELEMENT  CJDt 

m  —  —  —  —  —  —  —  —  —  — ——  f  ri 

-  VERSION  1.0 

— -  * 

-WRITTEN  bY  P . A . CUMO ALL ,  MARCH  19oD,  FOR  o.S.  AkM 
-(EUROPEAN  RESEARCH  OFFICE)  AND  Der ErtSE  NUCLEAR  AGEnCY 
-UNDER  CONTRACT  DAJA37-79-C-US4R . 


SUBROUTINE 

-UPDATE 


UPDATE 
DOMAIN  idom 


SUBROUTINE 

-UPDATE 


upot 

contacts  in  DOMAIN  idom 


subroutine 

-common 


VAR(NPAR) 

routine  for  ivar  & 


H  V  A« 


SUBROUTINE  XYCuR(A,X,Y) 

-RETURNS  CORNER  COORDS  ( X ,  Y)  FUR  Coki.ER  ARRAY  a  U 


SUBROUTINE  XYFO(A,X, Y,XD, YD) 

-RETRIEVE  GRID-POINT  COORDS  s»  VELOCITIES 


SUBRQUTii.e  ZCS 

—  TO  Rtl'DRi.  uEwGTri  (Z)  A*<D  Unix  vtCluR  (C/S) 
-OF  UnE  (XI ,  L  l ) ,  U2,  il ) 


SUBROUTINE  ZSTRSC A , IZOn ) 

-STRAlN-HATtS,  STRESSES  b  REriCt  G.R. 


FORCES  FOR  A  ZOwt 
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Appendix  III-4  Subroutine  Calling  Map 


UEEC 


-START 

-MON 


-10 

-TIDY 

-MATCH 

-10 

-IMSG 

-PROP 

-FIX 

-DAMP 

-DELB 


-10 


-SEP 


-MATCH 


-NEXTCP 
-DELC 
-LOSE 
-DELLST 
CHANGE  -MATCH 


-NEXTD  -DELLST  -LOSE 


-FIND 


-PRINT  -MATCH 
-OUTSD 
-OUTFDZ 

-OUTFDG  -OUTG 
-OUTCP 
-NEXTCP 
-NEXTD 

-NEXTP  -NEXTCP 


-CREATE 


-FIND 
-MATCH 
-GEO 
-CRAD 

-TOLSET  -C20 


-INERT 


-SPLIT 


-FIND 

-MOVE 

-GEO 

-NEXTCP 

-NEXTP 

-MAKEP 

-MAKEC 

-MAKEB 

-CRAD 

-COORD 


-TOLSET 

-UPDATE 


-INERT 


-NEXTCP 


-NEXTP 

-CC 

-CE 


-GEO 

-UPDT 


-NEXTCP 

-ZCS 

-REL1 

-ZCS 

-REL2 


-JUMPB 

-NEXTD 

-NEXTP 

-NEXTCP 

-NEXTDS 

-NEXTCP 

-QUERY 

-JUMPB 

-NEXTD 

-NEXTP 

-NEXTCP 

-NEXTP 

-NEXTCP 

-CE 

-RE  LI 

-ZCS 

-REL2 

-REL2 

-NEXTCP 

-CC 

-ZCS 

-BRIO® 

-FIND 

-NEXTCP 

-COORD 

-NEXT? 

-NEXTC 

-CC 

-ZCS 

-CE 

-REL1 

-ZCS 

-REL2 

-PUSH 

-FIND 

-APLOT 


-INI 


-PLOTS 

-ARC 


-PLOTS 


-INERT 


-CRAD 
-<Z0 
-MAKEB 

-TOLSET  -CZO 
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-CYCLE 


-NOTION  -CVEL 
-PVEL 

-STRAIN  -XYCOR 

-STRESS  -CL1 
-CL  2 
-CL  3 
-CL  4 
-CL5 

-FDEF  -CP MOT 

-XYFD 
-PVELFD 
-NEXTCP 
-2CS 
-REL2 
-CVELFD 


-2STRS  -XYFD 

-FORD  -NEXTP  -NEXTCP 

-CC  -2CS 

-CE  -REL1 

-ZCS 
-REL2 

-DELC  -NEXTD 


-DELLST 

-LOSE 

-CLJ2 

-CLJ3 

-CLJ4 

-CLJ5 

-PUTFCE 

-PUTSTR 

-FDC  -PUTFDF 

-FDE  -PUTFDF 

-MVFOR  -NEXTP  -NEXTCP 
-NEXTCP 
-MVBAK  -NEXTCP 

-NEXTP  -NEXTCP 

-ACCUM 

-NEXTP  •NEXTCP 

-NEXTD 

-UPDATE 

-GEN  -FIND 

-MATCH 


Most  "functions"  and  the  commoner  subroutines  which  are  listed  below  have 
been  left  out  of  the  "calling  map": 


Function* 

Subroutines 

(33  TR 

PUTI 

ICALB 

PUTP 

IDOMP 

POTR 

1GET 

VAR 

1GETGP 

IVAR 

QBC 

QTEST 

RANGE 

RVAR 
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APPENDIX  IV  -  Guide  to  Program  Changes 


IV- 1  Parameters 

All  Fortran  variables  initialised  in  the  BLOCK  DATA  subroutine  FRED 
are  constant  during  a  UDEC  run,  and  the  DATA  statements  may  be  replaced  with 
PARAMETER  statements  if  these  are  supported  by  the  compiler.  This  should  lead 
to  an  increase  in  efficiency. 

IV- 2  Storage  of  Integers  and  Real  Numbers 

Program  UDEC,  as  written,  assumes  that  a  real  Fortran  variable  occupies 
the  space  of  two  integers.  This  correspondence  is  only  important  for  quantities 
stored  in  the  main  array  IA (  ) .  If  UDEC  is  to  be  run  on  a  machine  having  a  dif¬ 
ferent  convention  for  storing  variables,  the  offsets  defined  in  the  BLOCK  DATA 
subroutine  will  need  to  be  changed.  The  offset  of  a  variable,  referred  to  a 
particular  data  array,  is  defined  as  the  number  of  words  from  the  start  of  the 
data  array  to  the  location  of  the  variable,  where  a  "word"  cam  be  integer  or  real, 
depending  on  the  type  of  variable ,  The  numbering  for  integers  starts  at  0,  and 
the  numbering  for  reals  starts  at  1.  As  am  exaunple,  the  offsets  are  given  below 
for  the  corner  data  array  for  the  standard  program  UDEC,  and  for  the  program  as 
it  would  be  set  up  on  a  machine  in  which  integers  occupy  the  same  space  as  reals. 
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Variable 


For 

2  integers 
s  1  real 


For 

1  integer 
=  1  real 


INTEGER  REAL 
COUNT  COUNT 


I _ I 

INTEGER  REAL 
COUNT  COUNT 


j  <head  code> 

(integer) 

0 

(1) 

0 

(1) 

|  KL 

(integer 

1 

1 

(2) 

?  KR 

II 

2 

(2) 

2 

(3) 

;  KNB 

II 

3 

3 

(4) 

‘  KXP 

(real) 

(4) 

3 

(4) 

5 

KYP 

It 

(6) 

4 

(5) 

6 

|  KXCP 

II 

(8) 

5 

(6) 

7 

|  KYCP 

II 

(10) 

6 

(7) 

8 

j  KRAD 

It 

(12) 

7 

(8) 

9 

KXDP 

II 

(14) 

8 

(9) 

lo 

KYDP 

II 

(16) 

9 

(10) 

11 

KGP 

(integer) 

(18) 

(10) 

11 

(12) 

numbers  not  in  parentheses  are 
offsets  used  in  program;  other  numbers 
are  included  to  show  the  full  sequence. 

IV- 3  Non-standard  I/O  and  other  Operations 

All  non-standard  input/output  is  done  in  subroutine  10,  and  consists 
of  OPEN  and  CLOSE  calls  for  the  input  and  output  files.  These  calls  may  be 
replaced  by  equivalents  if  another  computer  is  to  be  used. 

The  DECODE  statements  in  subroutines  IVAR  and  RVAR  may  be  non-standard 
on  some  machines,  and  should  be  replaced  as  necessary.  No  other  ENCODE  or 
DECODE  statements  are  used  in  the  program. 


data  array 
for  comers 

(see  Appendix  II I -2) 
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Subroutine  VAR  packs  characters  into  array  IBUF,  which  is  declared  a 

BYTE  array.  On  some  machines  this  operation  may  have  to  be  done  with  an  ENCODE 

* 

statement,  if  BYTE  or  INTEGER  1  variables  are  not  allowed. 

* 

All  INTEGER  2  statements  may  be  omitted  if  they  refer  to  single 
variables;  statements  that  refer  to  arrays  should  be  replaced  by  DIMENSION 
statements . 

The  statement  INCLUDE  ' COMMON. FTN'  inserts  at  that  point  in  the  program 
the  parameter  blocks  and  main  common  block.  Many  computers  have  a  similar 
facility,  but  is  invoked  differently. 

IV- 4  User-supplied  Constitutive  Subroutines 

Dummy  subroutines  CL2 ,  CL3 ,  CL 4  and  CL5  may  be  replaced  by  real  sub¬ 
routines:  these  will  be  called  for  materials  with  constitutive  numbers  2,3,4, 
and  5  respectively.  Input  and  output  variables  sire  passed  in  common  block 
/ CLCOM/ ,  where  the  names  of  the  variables  have  the  following  meaning: 


INPUT 


components  of 
current  stress  tensor 


components  of 

strain  increment  tensor 


OUTPUT 


components  of  stress 
increment  tensor 


